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Introduction and Context 

There has been a recent movement towards the use of advanced statistical analysis to track 
and project player and overall team performance in many professional sports. In particular, the 
paper "An Economic Evaluation of the Moneyball Hypothesis" by Hakes and Sauer [1] provides 
an analysis of Major League Baseball (MLB) player performance and identified which individual 
player statistics were more crucial to team success. Subsequently Hakes and Sauer determined 
which player attributes were undervalued within the MLB market based on salary. 

We wish to conduct a similar analysis with data from the National Basketball Associaton 
(NBA). Our goal is to identify which attributes contribute more to overall team success and 
attempt to predict team success using winning percentage. In this report, we construct a mul- 
tiple linear regression model on 30 years of NBA team data to predict this winning percentage. 
After selecting the best subset of predictors, we diagnosed our model and data for the pres- 
ence of multicolinearity, influential observations, goodness of fit, and autocorrelation. Our final 
model suggests that defensive tactics like creating turnovers and reducing opponent's field goal 
percentage are the most important contributors to team performance. 

Data Collection 

We obtained all data from Basketball-Reference.com. Initially, we began with individual player 
data [2] , but we soon realized that we would have difficulty associating individual player statis- 
tics with overall team success given the large number of players belonging to each team. To 
get around this problem, we collected team statistics (presented in basketball-reference on a 
per season basis) with the intention of regressing on winning percentage [3]. For each team, 
we combined three datasets from basketball-reference: "Team Stats", "Opponent Stats", and 
"Miscellaneous Stats". We first collected a raw team data set from a single year of observa- 
tions. The size of the data set was problematic, as we had 66 potential predictors of winning 
percentage with only 30 observations. We needed to collect more data to increase our sample 
size, but we did not want to introduce potential correlations between the same NBA teams 
from consecutive seasons. So, we decided to pull data from seasons a decade apart. This would 
minimize any potential correlations between teams that have the same players on their roster 
for more than one season. Our final data set contains 109 data points. 

Data Trimming 

Before we began to analyze the dataset for outliers and correlated predictors, we trimmed the 
dataset of numerous predictors. First, we removed pre-calculated statistics from out dataset. 
Basketball-reference.com provides some statistics such as Pythagorean Wins and Simple Rating 
System, which are not observed data and are calculations meant to predict the exact response 
we were trying to regress on (Wins or Winning Percentage) . Additionally, we trimmed out data 
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that was uniform across the entire dataset and all observations such as number of games played. 
This is because these observations do not add any additional predictive power as they are the 
same across all teams. 

There are multiple observations within the dataset that are likely to be highly correlated 
because they are interrelated, such as Free Throws Made (FT), Free Throws Attempted, Free 
Throw Rate, etc. However, we chose not to preemptively remove these variables. We did not 
know which variables would be the best at predicting our response, so we decided to remove 
them from the model later (see our Multicolinearity section). 



Model Selection 
Initial Selection 

We began by running a stepwise procedure with both AIC and BIC, regressing on Winning 
Percentage (Win. Pet.). The AIC stepwise procedure produced a model with an AIC of 661.362, 
an R 2 value of 0.936, and 31 predictors, while the BIC stepwise procedure produced a model 
with a BIC of 735.927, an R 2 of 0.927, and 25 predictors (26 estimated f3 — coef ficients). The 
discrepancy between the number of predictors makes sense as the penalty term for BIC makes 
that model selection criteria favor smaller models with fewer predictors. Since the R 2 values 
are almost identical, we chose to use the model produced by BIC selection criteria as it is more 
parsimonious. 



Multicollinearity 

With such a large number of predictors (and with many of them related), we knew there could 
be potential multicollinearity problems with our model fit. In the heat map of the correlation 
matrix associated with predictors in our BIC model (Figure 2), we saw numerous dark-colored 
cells indicating a large amount of correlation between many predictors. The lower left corner 
in particular indicates that those variables had large correlation coefficients. 
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Figure 1: Heat Map of Stepwise BIC Model 



Figure 2: Varience Inflation 
Factors of each predictor 



This was an indication that there was potential multicollinearity within our model. After 
calculating the Variance Inflation Factor for every predictor within our model, we found many 
predictors with VIFs above the standard threshold of 10 (figure 2, bolded). Since we had many 
predictors within our model (many of them measuring similar observations), we decided to run 
an iterative procedure where: 

1. The predictor with the largest VIF value was removed from the model 

2. VIF for each predictor was calculated for new model 

3. Process was continued until all VIF values were below 10 



2 



In the end, we removed seven predictors from our model: FG, MOV, FT, Opp. FG%, PTS.G, 
FTr, and Opp.TOV. We believe this was the correct decision because many of the predictors 
remaining in the model still accounted for the removed predictors (Opp.TOV% remained in the 
model, which is directly related to OPP.TOV, FG% remained in the model, which is direclty 
related to FG, etc.). We decided on this procedure instead of using Principle Component 
Analysis because we wanted to mantain the interpretive power of our model. 

Our final model (Figure 3) contains 18 predictor variables, and has an R 2 value of 0.8361. 
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Figure 3: Model after Removing Correlated Predictors 



The corresponding heat map for this model (figure 4) is an indication that the correlation 
issues have been resolved, with no dark-colored cells other than the diagonal. Additionally, 
none one of the associated varience inflation factors were greater than 10 (figure 5). 
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Figure 5: Variance Inflation 
Figure 4: Heat Map without Correlated Predictors Factor of each predictor 

Model Diagnostics 
Initial Examination 

The normal probability plot (Figure 6) is roughly linear. However, the sample quantiles deviate 
from their theoretical values at both the lower and upper ends of the plot. This indicates that 
some normality assumptions may have been violated along with a higher than average variance. 
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Normal Q-Q Plot 




-2-1012 
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Figure 6: Normal Probability Plot 

The residual plot (Figure 7) reveals some problems with our model fit. Though the residuals 
are fairly uniformly distributed in the horizontal band around zero, there are quite a few very 
large and very small residuals, indicating that there are potential outliers in our data set. These 
outliers give the residuals a slight funnel shape. We decided to conduct tests to see if there are 
any particularly influential observations. 



Model Residual Plot 
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Figure 7: Residual Plot 

Outlier Pruning 

We checked for influential observations using DFFITS as a measure. Our dataset is small 
with 109 observations, so our threshold for influential observations was if the magnitude of the 
DFFITS value is greater than one. 

After calculating DFFITS, we found 12 observations that were deemed influential. Since all 
12 observations came from the same season, which was an era where the style of play was much 
different than other eras, we decided that it would be a good idea to remove these outlying 
observations from our data set. We sill have 97 observations within our dataset, which is 
sufficiently large. 
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Estimate Std. Error t value Pr(>|t|) 
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Figure 8: Model with Outlying Observations Removed 



Predictor Reselection 

Note that many of the predictors in our model are insignificant when taken individually (Figure 
8) after we have removed the outlying observations. The t-tests for ORB, Opp. 3PA, Opp. 2p., 
TOV., FT.FGA, Opp.DRB., and Opp.FT.FGA all yielded p- values above a = 0.05. To test 
whether these predictors belong in the model when taken together, we performed the following 
partial F-test: 



Hq : PoRB = fiopp.ZPA = Popp.2p. = PtOV. = ^FT.FGA = ^Opp.DRB. = Popp. FT.FGA = 0 

Ha ■ At least one (3 coef ficient ^ 0 
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Under the null hypothesis, F* should follow F(7, 87). The p- value for our observed statistic 
(1.18) is 0.32, implying that we do not have sufficient evidence to reject Hq and that these 
predictors should be removed. 

Once we refit the final model with the new data, we get a much better normal probability 
plot that is almost completely linear, as well as much better residual plot with no patterning or 
outlying residuals (Figures 9 and 10). Additionally, our model fit has improved in the absence 
of outliers and insignificant predictors: our new R 2 is 0.936. 
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Figure 9: Normal Probability Plot Figure 10: Residual Plot 
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Transformation 



The diagnostic plots indicate that now the model is appropriate for the data, but we also ran a 
Box-Cox transformation just to ensure that no further improvements can be made to the model 
(Figure 11). 
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Figure 11: Box-Cox Plot 
As shown by the Box-Cox plot, no transformations are needed to improve the model fit. 

Autocorrelation 

Because our team data comes from successive NBA seasons, there is the possibility that corre- 
lations exist through time in our predictors. We attempted to combat this problem by spacing 
out each generation of teams by 10 years. When win percentage is plotted against teams (sorted 
by season year), there are no discernable patterns or trends in the data (figure 12). 



Win Percent per Team 
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Figure 12: Plot of Win. Pet Over Teams, sorted by Time 

Additionally, the autocorrelation function (ACF) plot indicates that there is little autocor- 
relation between our predictors (figure 13). The only lags with significant autocorrelation are 
1, 12, and 13. However, the magnitide of the ACF for these lags is barely beyond the 95% 
confidence threshold for autocorrelation, implying that the relationships are tenuous at best. 
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Series Win.Pct. 
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Figure 13: ACF Plot 

To test whether the true autocorrelation is greator than 0, we performed the following 
Durbin- Watson test with a significance threshold of a = 0.01: 

Hq : autocorrelation (p) = 0 
Ha '■ autocorrelation (p) > 0 

Our observed test statistic was D = 1.73, which translates to a p- value of 0.059. Because 
this is larger than our a, we fail to reject Hq and conclude that the true autocorrelation in our 
data is not different from 0. 

Results and Conclusion 

The Final Model 

Our final model is: 

WiruPct. = - 111.03 + 185.2x FG . - 0.02x DRB - 0.05x TO v ~ 0.02x Opp .fta 

+ O.Olxopp.PF + 307.99xt5. + 1.35x 0 rb. ~ 343.34xo PP . e FG. + 4.39xo PP .TOV. 
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Figure 14: Our Final Model with 9 Predictors 
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Model Verification 

We performed the following global F-test to validate our entire model: 

Hq : All (3 coef ficients = 0 

Ha '■ At least one (3 coef ficient ^ 0 



Source Df SS MS F 

Regression 9 19807.74 2200.86 

Error 87 1350.4 15.52 

Total 96 21158.14 141.8 

Under the null hypothesis, F* should follow F(10,87). The p-value for F* = 141.8 is 
2.2 x 10~ 16 , giving us sufficient evidence to reject the null hypothesis and conclude that at least 
one /3-coefficient is nonzero. Thus, our model fit is appropriate. 

Discussion 

The coefficients in our final model had both expected results and a few surprises. Due to the 
scaling of the percentage variables, the shooting statistics were calculated between the range of 
0 and 1, while winning percentage and rebound percentages ranged from 0 to 100, so we have 
to be careful when interpreting the size of the effect for each predictor. 

With this in mind, the predictor variable with the highest coefficient was actually opponent 
turnover percentage. A one percent increase in xopp.tumover% was associated with a 4.4 percent- 
age point increase to predicted winning percentage, on average. This makes sense as it would 
give a team more possession of the ball, having more chances to score while the opponent has 
less. Opponent's effective field goal percentage had the second largest coefficient. A one per- 
cent increase in opponent effective field goal percentage, which, due to the scaling, corresponds 
with a .01 increase in the predictor variable, is predicted to decrease winning percentage by 
approximately 3.43 percent. 

The most signifcant predictors affecting team performance concerned defensive statistics, 
such as reducing an opponent's effective field goal percentage and increasing their turnover rate. 
Thus a team's performance on defense is more indicative of their success. Using this model, in 
order to maximize winning percentage a team must focus on two aspects of basketball: 

1. Creating turnovers in order to maximize possession and thus shooting opportunities. 

2. Reducing the efficiency of opponent's shots 

These results seem to confirm the common basketbal adage: "Defense wins championships". 
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R Code 



#Load and attach data (Large Team Data Trimmed. csv) 
data = read. csv(file. choose () , header = TRUE) 
attach(data) 

#Model selection 

#Fit model with all explanatory variables included 

lm . f ull = lm (Win . Pet . "MP+FG+FGA+FG . +X3P+X3PA+X3P . +X2P+X2PA+X2P . +FT+FTA+FT . +ORB+DRB+TRB+AST+STL+BLK+TOV+PF+PTS+PTS . G+Opp . MP 
+0pp . FG+Opp . FGA+Opp . FG . +0pp . 3P+0pp . 3PA+0pp . 3P . +0pp . 2P+0pp . 2P . +0pp . FT+Opp . FTA+Opp . FT . +0pp . ORB+Opp . DRB+Opp . TRB+Opp . AST+Opp . STI 
+0pp . PTS+Opp . PTS . G+Age+M0V+S0S+0Rtg+DRtg+Pace+FTr+X3PAr+TS . +eFG . +T0V . +0RB . +FT . FGA+Opp . eFG . +0pp . TOV . +0pp . DRB . +0pp . FT . FGA) 

#A11 subset selection using stepwise procedure and AIC value 
step (lm. full, direction="both" , k=2) 

lm . best . aic=lm(f ormula=Win . Pet . "FG+FG . +X3P . +X2P . +FT+FTA+ORB+DRB+STL+TOV+PTS . G+Opp . FG+Opp . FGA+Opp . FG . +0pp . 3P+0pp . 3PA 
+0pp . 2P . +0pp . FTA+Opp . DRB+Opp . TOV+Opp . PF+MOV+FTr+TS . +T0V . +0RB . +FT . FGA+Opp . eFG . +0pp . TOV . +0pp . DRB . +0pp . FT . FGA) 

summary (lm.be st . aic) 

#A11 subset selection using stepwise procedure and BIC value 
step (lm. full, direction="both" , k=log(109)) 

lm . best .bic=lm(f ormula=Win . Pet . "FG+FG . +FT+FTA+ORB+DRB+STL+TOV+PTS . G+Opp . FG . +0pp . 3PA+0pp . 2P . +0pp . FTA+Opp . TOV+Opp . PF 
+M0V+FTr+TS . +T0V . +0RB . +FT . FGA+Opp . eFG . +0pp . TOV . +Opp . DRB . +0pp . FT . FGA) 
summary (lm.be st .bic) 



#Check correlation between variables in each of the models and display using heat map 
library (reshape2) 
library (ggplot2) 

variables=data . frame (MP , FG , FGA , FG . , X3P , X3PA , X3P . , X2P , X2PA , X2P . , FT , FTA , FT . , ORB , DRB , TRB , AST , STL , BLK , TOV , PF , PTS , PTS . G , 
Opp . MP , Opp . FG , Opp . FGA , Opp . FG . , Opp . 3P , Opp . 3PA , Opp . 3P . , Opp . 2P , Opp . 2P . , Opp . FT , Opp . FTA , Opp . FT . , Opp . ORB , Opp . DRB , Opp . TRB , Opp . AST , C 
Opp. TOV, Opp.PF, Opp. PTS, Opp. PTS. G, Age, MOV, SOS, ORtg.DRtg, Pace, FTr ,X3PAr ,TS . , eFG . , TOV. , ORB. , FT. FGA, Opp. eFG. , Opp. TOV. ,0pp. DRB. ,C 
varcor = cor (variables) 



Finalvars=dat a . frame (FG , FG . , FT , FTA , ORB , DRB , STL , TOV , PTS . G , Opp . FG . , Opp . 3PA , Opp . 2P . , Opp . FTA , Opp . TOV , Opp . PF , MOV , FTr , TS . , TOV . , ORI 
finalvarcor =cor (Finalvars) 



aic_var iables=dat a . frame (FG , FG . , X3P . , X2P . , FT , FTA , ORB , DRB , STL , TOV , PTS . G , Opp . FG , Opp . FGA , Opp . FG . , Opp . 3P , Opp . 3PA , 
Opp. 2P. , Opp. FTA, Opp. DRB, Opp. TOV, Opp.PF, MOV, FTr, TS. , TOV. , ORB. , FT. FGA, Opp. eFG. ,0pp. TOV. ,0pp. DRB. ,0pp. FT. FGA) 

aicvarcor = cor (aic_variables) 

#correlation matrix of full model 
regular_cor_melt = melt (varcor) 

qplot(x=Varl, y=Var2, data=regular_cor_melt , f ill=value,geom="tile")+ scale_f ill_gradient2(limits=c (-1 , 1)) 

#correlation matrix of best aic model 
aic_cor_melt=melt (aicvarcor) 

qplot (x=Varl , y=Var2, data=aic_cor_melt , f ill=value,geom="tile")+ scale_f ill_gradient2 (limits=c (-1 , 1)) 

(♦correlation matrix of best bic model 
cor_melt = melt (finalvarcor) 

qplot(x=Varl, y=Var2, data=cor_melt , f ill=value,geom="tile")+ scale_f ill_gradient2(limits=c (-1 , 1)) 



#We decided to go proceed with the bic model (discussion in paper) 



#Multicollinearity 

SEliminating correlated variables from the best BIC model 1 at a time untill all Variance Inflation factors are below 10 
Final varmatrix = cbind (FG , FG . , FT , FTA , ORB , DRB , STL , TOV , PTS . G , Opp . FG . , Opp . 3PA , Opp . 2P . , Opp . FTA , Opp . TOV , Opp . PF , MOV , FTr , TS . , TOV . 
, ORB. , FT. FGA, Opp. eFG. ,0pp. TOV. ,0pp. DRB. ,0pp. FT. FGA) 
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r. squared = rep(NA, 25) 
for (i in 1:25){ 

dat =cbind (data. frame (Y=Finalvarmatrix [ , i] ) , as . data. frame (Finalvarmatrix [ , -i] ) ) 
lm.fit = lm(Y ~ ., data=dat) 
r.squared[i] = summary(lm.fit)$r. squared 
} 

#Check previous Variance Inflation Ratior and reduce until all are below 10 
1 / (1-r. squared) 

# remove FG 

Finalvarmatrix = Finalvarmatrix [, -1] 
r. squared = rep(NA, 24) 
for (i in 1:24){ 

dat =cbind(data. f rame (Y=Finalvarmatrix [ , i] ) , as . data. frame (Finalvarmatrix [ , -i] ) ) 
lm.fit = lm(Y " ., data=dat) 
r.squared[i] = summary (lm . fit) $r . squared 
} 

# remove MOV 

Finalvarmatrix = Finalvarmatrix [, -15] 
r. squared = rep(NA, 23) 
for (i in 1:23){ 

dat =cbind(data. f rame (Y=Finalvarmatrix [ , i] ) , as . data. frame (Finalvarmatrix [ , -i] ) ) 
lm.fit = lm(Y " ., data=dat) 
r.squared[i] = summary (lm . fit) $r . squared 
} 

# remove FT 

Finalvarmatrix = Finalvarmatrix [, -2] 
r. squared = rep(NA, 22) 
for (i in 1:22){ 

dat =cbind(data. f rame (Y=Finalvarmatrix [ , i] ) , as . data. frame (Finalvarmatrix [ , -i] ) ) 
lm.fit = lm(Y " ., data=dat) 
r.squared[i] = summary(lm.fit)$r. squared 
} 

# remove 0pp. FG. 

Finalvarmatrix = Finalvarmatrix [,-8] 
r. squared = rep(NA, 21) 
for (i in 1:21){ 

dat =cbind(data. f rame (Y=Finalvarmatrix [ , i] ) , as . data. frame (Finalvarmatrix [ , -i] ) ) 
lm.fit = lm(Y " ., data=dat) 
r.squared[i] = summary(lm.fit)$r. squared 
} 

# remove PTS.G 

Finalvarmatrix = Finalvarmatrix [, -7] 
r. squared = rep(NA, 20) 
for (i in 1:20){ 

dat =cbind(data. f rame (Y=Finalvarmatrix [ , i] ) , as . data. frame (Finalvarmatrix [ , -i] ) ) 
lm.fit = lm(Y " ., data=dat) 
r.squared[i] = summary(lm.fit)$r. squared 
} 

#remove FTr 

Finalvarmatrix = Finalvarmatrix [, -12] 
r. squared = rep(NA, 19) 
for (i in 1:19){ 

dat =cbind(data. f rame (Y=Finalvarmatrix [ , i] ) , as . data. frame (Finalvarmatrix [ , -i] ) ) 
lm.fit = lm(Y ~ ., data=dat) 
r.squared[i] = summarydm.f it)$r . squared 
} 

# remove 0pp. TOV 

Finalvarmatrix = Finalvarmatrix [,-10] 
r. squared = rep(NA, 18) 
for (i in 1:18){ 

dat =cbind(data. f rame (Y=Finalvarmatrix [ , i] ) , as . data. frame (Finalvarmatrix [ , -i] ) ) 
lm.fit = lm(Y " ., data=dat) 
r.squared[i] = summarydm.f it)$r. squared 
} 

# get the vif to verify that all are below 10 
1/ (1-r. squared) 



# redo our model using these uncorrelated predictors 

new . best . bic . model=lm (f ormula=Win . Pet . "FG . +FTA+0RB+DRB+STL+T0V+0pp . 3PA+0pp . 2P . +0pp . FTA+Opp . PF+TS . +T0V . 
+0RB . +FT . FGA+Opp . eFG . +0pp . TOV . +0pp . DRB . +0pp . FT . FGA) 
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new.vars = cbind(FG. , FTA , ORB , DRB , STL,T0V,0pp.3PA,0pp.2P,0pp.FTA,0pp.PF,TS. ,T0V. ,0RB. , FT .FGA,0pp . eFG . 
,0pp. TOV. ,0pp. DRB. , Opp.FT.FGA) 

#Display correlation matrix of newest model using a heat map 

qplot (x=Varl , y=Var2, main='Model Correlation Matrix Heatmap', data=melt (cor (new. vars) ) , fill=value, 
geom="tile")+ scale_f ill_gradient2(limits=c (-1 , 1)) 



#Transf ormations 

#Determine if response variable can be predicted better after transformation 
library (MASS) 

boxcox (new . best . bic . model) 
#0utliers 

# initial diagnostic plots indicate that we have some outliers 
qqnorm(resid (new .best .bic .model) ) 

qqline(resid(new. best .bic .model) ) 

plot (r standard (new. best .bic .model) "Win. Pet . ,ylab="Studentized Residuals", xlab="Win Percent", main="Model Residual Plot") 

# check for outliers using DFFITS. Our data is medium/small so a critical value of 1 was used 
dfit = abs (df fits (new. best .bic .model) ) 

check <- c() 
ford in 1:109){ 

if else(df it [i] > 1, check<-c(check, i), NA) 
} 

check 

#remove these outliers 

to. remove = c("3", "6", "11", "13", "14" , "15" , "21" , "23" , "24" , "25", "27", "30") 
data. new = data[! rownames (data) '/.in'/, to. remove,] 
attach (data. new) 

#fit data with removed outliers 

new . best . bic . model=lm (f ormula=Win . Pet . ~FG . +FTA+0RB+DRB+STL+T0V+0pp . 3PA+0pp . 2P . +0pp . FTA+Opp . PF+TS . + 
TOV.+ORB.+FT.FGA+Opp.eFG.+Opp.TOV.+Opp.DRB.+Opp.FT.FGA, data = data. new) 

#Verify that residual plot does not have any noticeable patterns/features 

plot (r standard (new. best. bic. model)~Win. Pet. ,ylab="Studentized Residuals", xlab="Win Percent", main="Model Residual Plot") 



#Autocorr elation 

#Check for autocorrelation using acf method 
library (zoo) 
library (lmt est) 

plot (Win. Pet. , xlab = "Team", ylab="Win Percent", main="Win Percent per Team", type="b") 
acf (Win. Pet.) 

dwtest (new. best .bic .model) 

#Global F test 

anova (new. best .bic .model) 

atbl = anova(new .best .bic .model) 
atbl[,3] 

ssr = atbl[,3] 
ssr = ssr [-19] 
ssr = sum(ssr) 
msr = ssr/18 
mse = abtl[19,3] 
Fstar = msr/mse 

#F-statistic: 72.79 on 18 and 78 DF, p-value: < 2.2e-16 
#Reject null hypothesis 

#At least one BetaCoeff does not equal 0 
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#Partial F-test for insignificant variables 

reduced. best. bic. model = lm(formula = Win. Pet. " FG.+ DRB + TOV + Opp. FTA + Opp. PF + TS . + ORB. + 
Opp. eFG. + Opp. TOV., data = data. new) 

anova (reduced. best .bic .model , new .best .bic .model) 

#Null hypothesis not rejected- drop variables from model 

#diagnosis of new model 
qqnorm(resid(reduced.best .bic .model) ) 
qqline(resid(reduced.best .bic .model) ) 
boxcox (reduced . best . bic . model ) 

plot(rstandard(reduced. best. bic. model)~Win. Pet. ,ylab="Studentized Residuals" , xlab="Win Percent" , 
main="Model Residual Plot") 

anova (new. best .bic .model) 

Ftbl = anova (reduced. best .bic .model) 
Ftbl[,3] 

ssr = Ftbl[,3] 
ssr = ssr [-10] 
ssr = sum(ssr) 
msr = ssr/9 
mse = Ftbl [10,3] 
Fstar = msr/mse 

#F-statistic: 141.8 on 9 and 87 DF, p-value: < 2.2e-16 
#Reject null hypothesis 

#At least one BetaCoeff does not equal 0 
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